#include <stdio.h>
#include "./fasp/rtc.h"
#include "./fasp/vector.h"
#include "./fasp/qmatrix.h"
#include "./fasp/msolver.h"

QMatrix *Aum, *Avm, *Awm, *Apm;
Vector *bu, *bv, *bw,*bp;
Vector *xu, *xv, *xw,*xp;

void mem_alloca_u_(int *nbelements)
{
    Aum = malloc(sizeof(QMatrix));
    bu = malloc(sizeof(Vector));
    xu = malloc(sizeof(Vector));
    Avm = malloc(sizeof(QMatrix));
    bv = malloc(sizeof(Vector));
    xv = malloc(sizeof(Vector));
    Awm = malloc(sizeof(QMatrix));
    bw = malloc(sizeof(Vector));
    xw = malloc(sizeof(Vector));

    Q_Constr (Aum, "Momentum matrix x-component", *nbelements, False, Rowws, Normal, True);
    Q_Constr (Avm, "Momentum matrix y-component", *nbelements, False, Rowws, Normal, True);
    Q_Constr (Awm, "Momentum matrix z-component", *nbelements, False, Rowws, Normal, True);
    V_Constr (bu, "Momentum source x-component", *nbelements, Normal, True);
    V_Constr (bv, "Momentum source y-component", *nbelements, Normal, True);
    V_Constr (bw, "Momentum source z-component", *nbelements, Normal, True);

    V_Constr (xu, "Velocity x-component at cell center", *nbelements, Normal, True);
    V_Constr (xv, "Velocity y-component at cell center", *nbelements, Normal, True);
    V_Constr (xw, "Velocity z-component at cell center", *nbelements, Normal, True);
}


void mem_alloca_p_(int *nbelements)
{
    Apm = malloc(sizeof(QMatrix));
    bp = malloc(sizeof(Vector));
    xp = malloc(sizeof(Vector));
    Q_Constr (Apm, "Cont matrix p-component", *nbelements, False, Rowws, Normal, True);
    V_Constr (bp, "Cont source p-component", *nbelements, Normal, True);
    V_Constr (xp, "presure component at cell center", *nbelements, Normal, True);
}


//void mem_alloca(int *nbelements,double* apu,double *api,double *apn,double *bu)

void q_set_u_(int *element, int *numb, double *apu, double *apv, double *apw, double *bpu, double *bpv, double *bpw, int *api, double *apn_u,double *apn_v,double *apn_w)
{
   int i,j;

   Q_SetLen (Aum, *element, *numb + 1);
   Q_SetLen (Avm, *element, *numb + 1);
   Q_SetLen (Awm, *element, *numb + 1);

   Q_SetEntry (Aum, *element, 0, *element, *apu);
   Q_SetEntry (Avm, *element, 0, *element, *apv);
   Q_SetEntry (Awm, *element, 0, *element, *apw);

   for (j = 1; j <= *numb; j++)
   {
      Q_SetEntry (Aum, *element, j, api[j-1], apn_u[j-1]);
      Q_SetEntry (Avm, *element, j, api[j-1], apn_v[j-1]);
      Q_SetEntry (Awm, *element, j, api[j-1], apn_w[j-1]);
   }


   bu->Cmp[*element] = *bpu;
   bv->Cmp[*element] = *bpv;
   bw->Cmp[*element] = *bpw;

}

void q_set_p_(int *element, int *numb, double *app,double *bpp,int *api, double *apn){
   int i,j;

   Q_SetLen (Apm, *element, *numb + 1);
   Q_SetEntry (Apm, *element, 0, *element, *app);

   for (j = 1; j <= *numb; j++)
   {
      Q_SetEntry (Apm, *element, j, api[j-1], apn[j-1]);
   }
   bp->Cmp[*element] = *bpp;

}

void setuvw_(double *u, double *v, double *w, int *eleId){
    xu->Cmp[*eleId]=*u;
    xv->Cmp[*eleId]=*v;
    xw->Cmp[*eleId]=*w;
}

void solvematrix_fasp_u_(){

   double mtol=1.0E-8;
   int miter=10;
   double mres;
   int msolver=8;
   int mprecond=4;
   double mtime;
   int i;


  	SetRTCAccuracy (mtol);

        SolveMatrix (Aum, xu, bu, &miter, &mres, &mtime, msolver, mprecond, miter);
   	printf("\nMatrix U Number of iterations: %d Residual: %+E Time: %+E\n",
			miter, mres, mtime);

        SolveMatrix (Avm, xv, bv, &miter, &mres, &mtime, msolver, mprecond, miter);
   	printf("\nMatrix V Number of iterations: %d Residual: %+E Time: %+E\n",
			miter, mres, mtime);

        SolveMatrix (Awm, xw, bw, &miter, &mres, &mtime, msolver, mprecond, miter);
   	printf("\nMatrix W Number of iterations: %d Residual: %+E Time: %+E\n",
			miter, mres, mtime);
}

void output_uvw_(int *step){

	WriteMatrix (Aum, -1,1,*step);
	WriteVector (bu,1,*step);
	
	WriteMatrix (Avm, -1,2,*step);
	WriteVector (bv,2,*step);

	WriteMatrix (Awm, -1,3,*step);
	WriteVector (bw,3,*step);
}

void getuvw_(double *u, double *v, double *w, int *eleId){
    *u=xu->Cmp[*eleId];
    *v=xv->Cmp[*eleId];
    *w=xw->Cmp[*eleId];
}

void getpc_(double *p,int *eleId)
{
  *p=xp->Cmp[*eleId];
}

void solvematrix_fasp_p_(){

   double mtol=1.0E-9;
   int miter=100;
   double mres;
   int msolver=9;
   int mprecond=0;
   double mtime;
   int i;

   SetRTCAccuracy (mtol);
   SolveMatrix (Apm, xp, bp, &miter, &mres, &mtime, msolver, mprecond, miter);
   	printf("\nMatrix P Number of iterations: %d Residual: %+E Time: %+E\n",
			miter, mres, mtime);
}

void output_p_(int *step){

	WriteMatrix (Apm, -1,4,*step);
	WriteVector (bp,4,*step);
}


void mem_dealloca_u_()
{
   Q_Destr (Aum);
   Q_Destr (Avm);
   Q_Destr (Awm);
   V_Destr (bu);
   V_Destr (bv);
   V_Destr (bw);
}


void mem_dealloca_p_()
{
   Q_Destr (Apm);
   V_Destr (bp);
   V_Destr (xp);
}

